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On the basis of recent investigations, a newly developed analytical procedure is used for con- 
structing a wide class of localized solutions of the controlled three-dimensional (3D) Gross-Pitaevskii 
equation (GPE) that governs the dynamics of Bose-Einstein condensates (BECs). The controlled 3D 
GPE is decomposed into a two-dimensional (2D) linear Schrodinger equation and a one-dimensional 
(ID) nonlinear Schrodinger equation, constrained by a variational condition for the controlling po- 
tential. Then, the above class of localized solutions are constructed as the product of the solutions of 
the transverse and longitudinal equations. On the basis of these exact 3D analytical solutions, a sta- 
bility analysis is carried out, focusing our attention on the physical conditions for having collapsing 
or non-collapsing solutions. 

PACS numbers: 02.30.Yy, 03.75.Lm, 67.85.Hj 



I. INTRODUCTION 

About 85 years ago, the seminal work of Bose [T] opened up the study of the statistical properties of bosons in 
ultra-cold quantum systems. Bose's idea was further developed by Einstein [2], leading to the theoretical prediction 
of the condensation of atoms in the lowest quantum state below a certain temperature. The idea of Bose-Einstein of 
atom condensation in the ground state has been experimentally verified in a dilute gas composed of atoms [31 HI [51 E] ■ 
The dynamics of nonlinearly interacting bosons in ultra-low temperature gases is governed by the Gross-Pitaevskii 
equation (GPE), which is an extension of the nonlinear Schrodinger equation (NLSE) by including the confining 
potential and inter-atomic interactions [TJ. The GPE, without the external potential, admits localized solutions in the 
form of one-dimensional dark and bright solitons, as well as radially symmetric vortex structures. Nonlinear localized 
excitations involving BECs arise due to the balance between the spatial dispersion of matter waves and nonlincaritics 
caused by repulsive or attractive inter-atomic interactions in BECs. Recent observations [5J O E3 EQ HH E3 Ell E3 
conclusively demonstrated the existence of bright E3 EH EH ED and dark/grey [HI M EH EQ matter wave solitons 
and quantum vortices |16j . 

Although the area of investigations of localized solutions of the GPE is quite fascinating, most of the theoretical 
results deals with approximate solutions in 3D or in reduced geometries |17j . They are well supported by suitable 
numerical evaluations [18 and adequately compared with a very broad spectrum of experimental observations (for 
a review, see Ref.s [HI BO].) Nevertheless, this testifies that finding exact localized solutions of the 3D GPE in a 
trapping external potential well, and preserving their stability for a long time, is still a challenging task. In particular, 
one encounters serious difficulties in attempting to find soliton solutions in one or more spatial dimensions, although 
several kinds of solitons have been found using certain approximations [21 , 22]. This leads us to arrive at the conclusion 
that, in order to have exact soliton solutions in BECs, some sort of "control of the system" may be necessary. 

The very large body of experience suggests that interacting bosons constitute a nonlinear and nonautonomous 
system [23], for which coherent stationary structures (i.e. solutions of the 3D GPE) exist only if suitable time- 
dependent external potentials arc taken "ad hoc" ,24J. Therefore, the correct analysis of the system should include 
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a 'controlling potential' term in the GPE, to be determined self-consistently with the desired solutions (e.g. the 
localized solutions). This procedure may be, in principle, extended to an arbitrary 'controlled solution' with the 
appropriate choice of the controlling potential [25] . The controlling potential method (CPM) has been proposed in 
the literature, and used to find the multi-dimensional controlled localized solutions of the GPE [25 . Preliminary 
investigations |26j have suggested that control operations introduced by this method ensures the stability of coherent 
solutions against relatively small errors in experimental realizations of the prescribed controlling potential. This idea 
could be realized by techniques that involve lithographically designed circuit patterns, providing the electromagnetic 
guides and microtraps for ultracold systems of atoms in BEC experiments |27j . and by the optically induced 'exotic' 
potentials [25] . 

Another important aspect of solitons in BECs is the phenomenon of collective collapse/explosion, that has been 
predicted [22] and observed experimentally [251130]. In particular, this phenomenon seems to be dependent on the 
parameters of the BECs and on the confining or repulsive potential [22] . 

The stabilization and control of BECs in asymmetric traps have been investigated via time-dependent solutions 
of the GPE [31]. Stable condensates, with the limited number of 7 Li atoms with attractive interaction, have been 
observed in a magnetically trapped gas [32] . 

Recently, a mathematical investigation oriented towards the construction of 3D analytical solutions of the controlled 
GPE has been carried out [33] and applied to the construction of 3D exact localized solutions 34J. In Ref. [33 , 
it has been proven that, under the assumption of the separability of the external trapping potential well Vt ra p 
in the spatial coordinates [viz. Vt ra p{x,y, z,t) = V±(x, y, t) + V z (z,t), where V± and V z are referred to as the 
'transverse' and the 'longitudinal' potentials, respectively], and a suitable constrained variational condition for the 
controlling potential V con t r (i.e. the average over the 'transverse' x — y plane of V con t r is required to be a stationary 
functional of the BECs transverse profile), the factorized form of the solution of the 3D controlled GPE, in the form 
i)){x,y,z,t) = ip_i_(x,y,t)il} z (z,i), can be constructed, so that ip± and ip z satisfy a 2D linear Schrodinger equation and 
a ID nonlinear Schrodinger equation, respectively. 

In this paper, we apply the results of the above investigation [33] to develop a new analytical procedure for 
constructing a broad class of exact localized solutions of the controlled 3D GPE, with a parabolic external potential 
well. In particular, we extend our recent investigation [34 to a wider family of exact 3D localized solutions of the 
controlled GPE and perform an analysis that establishes the physical conditions and parameter regimes for having 
collapsing and non-collapsing localized solutions. In the next section, we formulate our problem and present the 
controlled GPE, and we briefly summarize the results found in Ref. |33j . In section III we apply these results to obtain 
localized solutions of the controlled GPE in the form of bright, dark and grey solitons for the longitudinal profile l^zl 2 , 
and in the form of the Hermite-Gauss functions for the transverse profile |?Aj-| 2 - We use the decomposition procedure 
of the controlled GPE suggested in Ref. [33] and solve the 2D transverse linear Schrodinger equation to obtain ip± 
in terms of Hermite-Gauss modes. Then, the ID longitudinal controlled NLSE is solved using a method based on 
the Madelung's fluid representation [33] [3B], which separates the NLSE into a pair of equations, composed of one 
continuity equation and one Korteweg-de Vries (KdV)-type equation. It is shown that the phase of the longitudinal 
wavefunction ip z has a parabolic dependence on the variable z. As the transverse and longitudinal equations are 
coupled both through the coefficient of the nonlinear term (in the longitudinal equation) and through the controlling 
potential, the consistency condition between the transverse and longitudinal solutions set up a relationship between 
the transverse and longitudinal restoring forces of the external trapping potential well. From the latter, the explicit 
spatio-temporal dependence of the controlling potential is self-consistently determined in each particular case. In 
section |IV| a detailed analysis of the dynamics of the above exact 3D solutions is developed both analytically and 
numerically. In particular, we study the properties of the controlled BEC states in the 2D case, showing that they 
feature breathing (oscillations of the amplitude and position) due to the oscillations of the perpendicular solution, 
as well as the oscillations in the parallel direction, arising from the initial displacement of the structure from the 
bottom of the potential well in the parallel direction. A stability analysis of controlled GPE structures is carried out 
in terms of a set of six parameters related through four equations. Therefore, two of them can be assumed arbitrarily. 
We discuss some examples of parameters corresponding to collapsing and non-collapsing 3D solutions. Finally, the 
conclusions are summarized in section [V] 



II. DECOMPOSITION OF THE CONTROLLED GROSS-PITAEVSKII EQUATION 

It is well known that the spatio-temporal evolution of an ultracold system of identical atoms forming a Bose Einstein 
condensate (BEC) in the presence of an external potential U ex t(v, t), within the mean field approximation, is governed 
by the 3D GPE, viz. 

lh l^ = "2^T y2 * + ^l*! 2 * + Uext{r,t) 9, (1) 
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where ^(r, t) is the wavefunction describing the BEC state, m a is the atom mass, Q is a coupling coefficient related to 
the short range scattering (s-wave) length a representing the interactions between atomic particles, i.e. Q = 4-7r?i 2 a/m a , 
and N is the number of atoms. Note that the short range scattering length can be cither positive or negative. Solitons 
have been observed in BECs of 7 Li atoms with a small scattering length a sa — 0.2nm, in correspondence of the 
following typical values of the parameters: N = 10 4 -10 5 at a temperature of 1-10 pK and a magnetic field ~ 400- 

600 g [El EI]. 

We assume that U ext is the sum of a 3D trapping potential well, say V tr ap, to confine the particles of a BEC, 
and a controlling potential, say V con t ri to be determined self-consistently. Furthermore, under this assumptions, we 
introduce the variable s = ct (c being the speed of light) and divide both sides of Eq.([l]) by m a c 2 , in such a way that 



Uext ^ 2 t) = V trap {v, s) + V contr (r, s) , (2) 



m a C 



and Eq. ([!]) can be cast into the form 

_ 

lh ~ 2 



dib \ 2 

i A c = -4r v ^ + [Vtra P (r, s) + V contr (r, s) + <#| 2 ] V , (3) 



where ip(r,s) = *&(r,t = s/c), \ c = Ti/m a c 2 is the Compton wavelength of the single atom of the BEC and q = 
NQ/mc 2 . 

In order to decompose Eq. ([3]), we briefly summarize the results of Ref. |33j . To this end, we first assume that 

V t rap(r,s) = Vx{r x ,s) + V z {z ) s), (4) 

where, in Cartesian coordinates, r = (x,y,z) and rj_ = (x,y) denotes, by definition, the 'transverse' part of the 
particle's vector position r and z the 'longitudinal' coordinate. Additionally, let us denote, in Cartesian coordinates, 
Vj^ = xd/dx + yd/dy. 

Let us suppose that -i/>^(r^, s) and ip z (z, s) are two complex functions satisfying the following 2D linear Schrodingcr 
equation 

i\c-£--H± > )l> ± = 0, (5) 



d 



•s 



where 



and the following ID NLSE 



where 



H± = -4vl+V x (r x ,s) !(.) 



H z = - > i^ + V z (z,s) + q 1D ( s )\^ z (^s)\ 2 + Vo, (8) 



2 dz 2 

respectively. In the latter, Vq is an arbitrary real constant and the function <7id(s) is defined as 



qiD(s) = q J d z r x ^J* ■ (9) 

Hereafter qm is referred to as the 'controlling parameter'. 

Furthermore, let us assume that the controlling potential depends in principle on rj^ through | | 2 , viz. Vcontr — 
Vcontr z,s), where px(rx,s) = \ipx(*x, s)\ 2 . 

Provided that 

V contr (rx,z,s) = [q 1D (s) - q\i;x(rx,s)\ 2 ] \ip z (z, s)\ 2 + V , (10) 
which makes stationary the functional (with respect to variation Spx of px', z and s play here a role of parameters) 

V[px; z, s] = J px{rx,s) V con tr {p±(r x, s) , z , s) d 2 rx , (11) 
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under suitable constraints provided by the normalization condition for ip±, viz. J pj_(rj_, s) dr± — 1, and by Eq. (j9j), 
where quj is thought as a given function of s, it can been shown |33j that the complex function 

iI>(t,s)=1>x(t x ,8)1> s (z,s) (12) 

is a 3D solution of the controlled Gross-Pitaevskii equation ^ . 

In summary, according to the results given in Ref. |33j . one may construct 3D solutions in the factorized form ( |12[ ) 
such that the 3D controlled GPE is decomposed into the set of coupled equations ^ and ^ plus the self-consistent 
expression of the controlling potential (10) coming from a constrained variational condition. Once Eqs. ^ and ^ 



are solved, ipx( r Xj s) and ip z (z, s) become known functions. Consequently, the explicit space and time dependence of 
V con tr is self-consistcntly determined, as well, i.e. the appropriate controlling potential corresponding to the controlled 



solution (12 1 is deduced. 

Note that, since ipx satisfies Eq. according to the definition (111, V represents the average of V cori tr in the 



transverse plane. The value of this average corresponds to the arbitrary constant Vo. Without loss of generality, we 
may put Vq = 0, viz. 

'<Pf f x1>lV eontr il> ± = 0. (13) 

In this way, among all possible choices of V con t r , we adopt the one which does not change the mean energy of the 
system (the average of the Hamiltonian operator in Eq. ^ is the same with or without V con t r ) and therefore minimizes 
the effects introduced by our control operation. 

In the next section, we will apply the results obtained here to the case of parabolic potentials, Vj_ and V z , to give 
exact 3D controlled localized solutions of Eq. ^ . 

III. EXACT LOCALIZED SOLUTIONS OF THE CONTROLLED 3D GPE WITH A 3D PARABOLIC 

POTENTIAL WELL 

Let us assume that Vj_(rj_, s) and V z (z, s) are the usual parabolic potential wells to confine the particle of a BEC, 
viz. 

Vx{rx,s) = \[ul{s)x 2 + ^( S )y 2 ] , (14) 

and 

V z (z,s) = ^ 2 z (s)z 2 , (15) 

where, in general, the frequency j — x,y,z, are supposed functions of time. The standard confining potential 
wells (along each direction) corresponds to the assumption that they are real quantity (positivity of their squares). 
However, our analysis can be extended to the case in which they are assumed imaginary (negativity of their squares). 
It follows that Eq. (|3| becomes 

iXc its = ~\ V ^ + \ ^ (S) X * + W * (S) ^ + W * {S) z2 ^ + q l^' 2 ^ + Vcontr y ' z ' s ^ ■ ( 16 ) 
According to the results and the assumptions of the previous sections, if we seek a solution in the factorized form 

i>{x, y, z, s) = ipx(x, y, s) ip z (z, s), (17) 



Eq. (16 1 can be decomposed into the following set of equations: 



1 Ac it + ¥ - \ [ul {s) x2 + ^ (s) y2] ^ = °' (18) 

. ^ d^^s) + | X_t0A _ 1 ^ {s) z 2^ g) _ qMs) ^ (Zj s)| 2 ^ ^ s) = Qi (lg) 



V contr (x,y,z,s) = [q 1D (s) — g|^j_(r^, s)| 2 ] \ip z (z,s)\ 2 . 



(20) 
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A. Solution of the transverse equation with a 2D parabolic potential well 



Equation (18 1 is readily solved, and its solution can be found in the standard literature, but we present it here for 
completeness. The general solution ip±(x,y,t) can be expressed as the superposition of Hcrmite-Gauss modes, viz. 

oo oo 

ip± (x, y, s) = ^2 ^2 an - 1 ^ x > n ( x ' s ) ^v,i (f i s ) ( 21 ) 
where a„j are arbitrary constants and 



n=o ;=o 



i> j>k (j,s) = [7r2 2k+1 (k\) 2 a 2 (s)\ 



exp 



ijj {s)j 2 
2A, 



4a 2 (s) ' 



Hi 



J 

\/2<t 3 {s) 



(22) 



where j — x,y. The perpendicular spatial scale <jj (i. e., the root of mean square) of the Hermite-Gauss functions 
( 22 1 satisfies the Ermakov-Pinney equation [STJ IMj 



d 2 crj (s) 
ds 2 



+ u 2 (s)a j (s) 



X 2 



4<r ( S ) 



= 0, 



and the phase functions jj(s) and <fij,k(s) are given by 

lj ( s ) 



o-j 0) ' 



<Aj,fc (s) = <t>j,kfi ~ -j ( 2k + !) 



ds 1 



(23) 

(24) 
(25) 



where 0j,fc,o is an arbitrary constant. 



B. Exact solution of the longitudinal NLSE with a ID parabolic trapping potential well 



In order to solve Eq. (19), we first observe that, according to Eqs. (f9j> , pl| and (22 1, the controlling parameter 
qm can be expressed in terms of the features of viz. 



,,_/M«),ff,W] 
(T^ (s) <J y (s) 



where T [(J x {s),a y {s) 1 s\ is a (relatively complicated) positive definite functional of o~ x and a y given by 

— 7= / — 7= y~] F n , m ,i,p[°x{s),ay(s)]H n (u)H m (u)Hi(v)H p (v 

with it = x/^2a x , v = yj\[2a y and 



n,m,l,p 



J r n,m,l,p[Ux{s),a y (s)] = 



y / 2 7l+m + i +p n\ ml Up] 



o »(0x,n,O-0x,m,o) 



o^'AM.i.O-^H.P-o) i 



(26) 



(27) 



(28) 



From the above equation, it is clear that, in general, qm depends in a non-trivial manner on a x (s) and cr y (s). 
However, in the simple case when the perpendicular solution contains only one Gauss- Hermite mode, say the (n, V)- 
mode, Eqs. (21 1 and (22 1, T becomes a real positive constant, viz., 



T[a x {s),a y {s)] = 
and, therefore, q\r> becomes 



H*(u)du f e- 2v Hf(v)dv _ 



/2 7T 



qiD(s) 



$n Si 



a x (s) o-y (s) 



V2n 
= <7iM s ) ■ 



S n Si — constant 



(29) 



(30) 
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1. Reduction of the ID NLSE to a KdV-like equation by means of the Madelung's fluid approach 



Equation (19 1 is a ID GPE with a time dependent parabolic potential. Its approximate solutions are well known in 



the literature, but we attempt here to find an exact solution which is also compatible with the ones of the transverse 
part, in such a way to give a solution of the full controlled 3D GPE (16 1. We seek ip z (z, s) using the following standard 
Madelung's fluid representation 



ip z (z, s) = y/p(z,s) 



exp 



iO(z, s) 

Xr 



which, after the substitution into Eq. ( 19 I and the separation of real and imaginary parts, yields 

2 



dO 
~ds~ 



dO 

dz 



tt< \ X " 1 d2 P U2 n 



and 



where 



dp 

ds 



d ( dO 



dz \^ dz 



= 0, 



(31) 



(32) 



(33) 



U(z,s) = 2 W z( s ) z 2 + qi D (s) p(z,s) 



(34) 



is the 'longitudinal potential energy' which is a functional of p. 

By differentiating Eq. ( |32| with respect to z and introducing the 'current velocity' V = dO/dz, a non-trivial series 
of transformations given in Ref.s [351 13S] allows us to obtain the following generalized KdV equation 



dV 



V 



dp 

ds 



co(s) - 



dV 

~ds 



dz 



dp 

dz 



dU 
dz 



2U 



dp\ 1 d 3 p 



dz 



Adz 3 



= 0, 



where co(s) is an arbitrary function of s. Making use of definition (34l, Eq. (35) becomes 
dV 



- P 



ds 



V 



dp 
ds 



co(s) 



dV 
ds 



dz 



d P .2/ \ 2/ \ 2 d P o„ d P , A c d3 P 



dz 



u z {s) zp-uj z (s) z" 



dz 



3qi D (s) p 



dz 4 dz 3 



(35) 



(36) 



Therefore, the system of equations (32 1 and (33) is now replaced by the system of equations (36 1 and (33) 



We look for a solution for p, assuming that V is a linear function of z, viz. 

V{z,s) = g{s)z + n{s) 
and this corresponds to seek a quadratic form of the solution for the phase 9(z, s), viz. 

6 (z, s) = O (s) + k (s) z + i g(s)z 2 , 



(37) 



(38) 



where the 'initial phase', 0o(s), the 'wavenumber', k(s), and the 'dispersive' term, g(s), are real functions of the 
time-like variable s. It is easy to see that the arbitrary function of s, co(s), appearing in Eq. (351, is proportional to 
Op (hereafter the prime stands for the first-order derivative with respect to s), i.e. cq(s) = — 0q(s). 

After substituting Eq. ( 37 1 into Eq. (36 1, the system of equations d33| and ( 36 1 can be cast into the form 



dp . .dp 

ir + {gz + K) ir +gp 

ds dz 



0. 



(39) 



and 



p{g' ' z + k') + 



OS 



-20d 



/ 2 

9 z 



2k' z 



dp 

dz 



2 2 2 dp dp Ac d 3 p 

UJ z Zp-u) z Z 3qiD p-K- + — ^7 

dz dz 4 dz A 



(40) 
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Then, substituting Eq. ( 39 1 into Eq. ( 40 1 we obtain 



dp A 2 d 3 p 



(41) 



(42) 



where 0' (s) = Q' (s) + k 2 (s)/2. To reduce Eq. (41 1 to the following KdV-like equation 

-20 o ( S )-- 3(?lc ( S )p- + T ^^O, 

we have to impose that the coefficients of (^p + z^j z and {^p + 2 zpfj are zero, namely we automatically find that 
g(s) satisfies the following Riccati's equation, viz. 



while k(s) is related with it through 
which is readily integrated as 



.g' + 5 2 + ^ 2 = 0, 
k(s) = k e-Zo'fW^, 



(43) 
(44) 
(45) 



where n is an arbitrary constant. 

We look for functions p(z, s) which satisfies simultaneously the KdV-like equation (42 1 and the continuity equation 
(39]). By using Eqs. (31 1,(p8|) and (43 1-(45 ), they allow us to construct also solutions of the longitudinal equation 
(196- 



To this end, under the coordinate transformation 

£ = £,{Z:S) = q 1D (s) z + R(s) 

s = s'(z, s) = s , 



(46) 



where R(s) is a real function, the system of Eqs. (|39| and (42 1 becomes 



Qid 



g)^-R) + R' + kq id 



dp dp 



and 



dp , a 2 2 ajV 



= 0, 



where the prime denotes differentiation with respect to s. 
To find solutions in the factorized form 



p(Z,s')=A(s')F(0, 



satisfying simultaneously (47 1 and (48 1, the following conditions have to be imposed 

Qid + 9 Qid = 

R' + nq 1D = 0. 



Consequently, Eqs. (47 1 and (48 1 become the following ordinary differential equations, respectively, 

A' + gA = 0, 

and 

dF 



29' 



dF A 2 2 d 3 F n 
-0. 



o _^ 3qiDAF _ + triD _ 



(47) 



(48) 



(49) 



(50) 



(51) 



(52) 
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We note that the first condition (50 1 and continuity equation (51 1 imply, respectively 

qiD(s') = q e-So'a(r)dr j 



and 



A(s') = A Q e-f«' g( - T)dT , 



(53) 



(54) 



where qo and Aq are integration constants, i.e. qo — qio(s' = 0), Aq = A(s / = 0). Additionally, by using solutions 
(53 1 and (45), the second condition (50 1 can be easily solved for R(s'), viz. 



R(s') = R - n qo / ds"e- 2 f<>" 9{T)d 



(55) 



where R = R(s' = 0). 



Furthermore, we also observe that, given the set of Eqs. (43 1, (451, (53 1 - (551, we can conveniently express the 



functions A(s), k(s) and u> z (s) in terms of the controlling parameter qio(s) (which is directly connected with the 
transverse part of the GPE solution), i.e. 



A(s) = — q 1D (s) , 

qo 



(56) 



k(s) = — qi D (s) , 

qo 



(57) 



R(s) = ~- / qi D (T)dr + R Q . 

qo Jo 



(58) 



<^(s) = -qiD(s) 



j2 r 



ds 2 



qiD(s) 



(59) 



In particular, Eq. (59) has been obtained by substituting the first of equations (|50| into Riccati's equation (43). It 



establishes a 'control condition' by the transverse part of the GPE solution on the longitudinal parabolic potential. 
In fact, it indicates which time dependence of uj z — lo z (s) has to be taken, provided that qm(s), given by Eq. (26 1, 
is solution of the first of Eqs. (50 1. Note that, according to Eq. (26), it results that 



So 



T [(Tx(s), (Ty(s)] s = Q 

a x (s = 0) a y (s = 0) 



(60) 



(Note that ^ r [cx{s),a y (s)] s _ does not coincides with !F[a x (s = 0),<r y (s = 0)].) To cast Eq. (52) as an equation with 
constant coefficients, we can choose the arbitrary function &' (s') proportional to e~ 2 f° 9 ^ dT and therefore 



(-),',(. s' I = %1 D { S ') 

qo 



(61) 



where 0q is an arbitrary constant, in such a way that, according to Eq. (451 and the definition of #o( s ') given above, 
we have 



Qp + K o/ 2 2 f J\ _ u 2 , ! 

2 q\D\ s ) — ~2 q\D\ s > ■ 
qo 



(62) 



with constant coefficients (stationary KdV equation 



dF 



By substituting Eqs. (|56j, ( |57| and ( |62[ ) in Eq. ( |52| , we finally obtain the following ordinary differential equations 

= 0. (63) 



dF \ 2 c q 2 d 3 F 



d£, 3 
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Let us now determine the phase Q(z, s). To this end, we first observe that Oo(s) can be obtained by integrating 
Eq. (f6lT), i.e. 



6' f s 

6o(s) = ¥'o + -2 / 9ii)( T ) dT : 
% Jo 



(64) 



where ipo is an integration constant. Without loss of generality, we can put: tpo = 0. Then, by using the first of 
conditions (50) and Eqs. (57) and (64 1, Eq. (38 1 can be easily cast into the form 



G(z,s) 



1 q'id(s) 



where 



2qi D (s) 

K o ql D (s) 



z - C(a) + A(s) , 



C(*) = 



io q'i D ( s ) 



(65) 



(66) 



(67) 



2. Soliton solutions 



As it is well known, Eq. (63 1 admits both localized and periodic solutions [321 HOI SI]- Typically, the latter are 
expressed in terms of Jacobian elliptic functions, whose suitable asymptotic limits of their parameters recover the 
localized solutions in the form of bright, dark and grey solitons. However, a very useful integration approach of Eq. 
( 63 1 that has a simple physical meaning of the solutions is the well-known Sagdeev's pseudo-potential method [HI [33] ■ 



(i). Bright solitons 



If F (and consequently ip z ) is a normalizable wave function, we can look for a bright soliton solution of Eq. (63 1 



which satisfies the following boundary conditions: F — > 0, for £ — > ±oo. This solution exists for 9 > , [36], i.e 

•2ft! ■ I \^ 



F(0 = f sech 2 I f— 1£ 



(68) 



Then, going back to the old variables, z and s and using Eqs. (491, (56l and (58l, we finally get the following bright 
soliton solution of Eq. ( 42 ) 



p(z, s) = —^-§-q\D{a) sech 2 
% 



A c |<7o 



■qiD(s) [z - z (s)} 



(69) 



where 



Zo(s) 



/ Qid{ t ) dr '- -Ho 
Jo 



(70) 



The positivity of p(z, s) implies that quj < and therefore q < and go < 0. On the other hand, if we require that 
ip z is normalized, i.e. 



J dz\ip z 



1. 



then 



1o _ I 2 \ • ?r [ (T ^( s )> cr a( s )] s= o 



Ox (s = 0) (s = 0) 



(71) 



(72) 
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It follows that, in the case of bright solitons, Eq. (67 1 becomes 



1 



ql D (r)dT + 



«o gig (a) 
?o i'id( s ) 



(73) 



(ii). Grey and dark solitons 



If F (and consequently ip z ) is a non-normalizable solution, we can look for dark or grey soliton solutions of Eq. 



(63 1, which satisfy the following boundary conditions: F — > Fa, for £ — > ±oo, where |i*o| < oo. These solutions are 

F(0=F [l-e 2 sech 2 (£/A )] , 



given by the general form 



where e is a real parameter, 



and 



F 



A = 



20' 



q A (e 2 - 3) 



\ 2 r, 2 

26' 



1 - 



We first observe that since Ao > 0, the following condition holds: 

0'o (e 2 - 3) > 0. 



Taking into account Eqs. (49) and ( 74 1-( 76 ) , we easily get: 



P(^s') 



28' Q q 1D (s') 
<Z 2 (e 2 -3) 



1 - e 2 sech 2 



26' 



1-3/e 2 A c |g | 



(74) 



(75) 



(76) 



(77) 



(78) 



which, going back to the variables z and s, can be cast into the form 

29' q 1D { S ) 



p{z,s) 



qli^-i) 



1 - e 2 sech 2 



26' q 1D (s) 



1-3/e 2 A C M 



(z-z (s)) 



(79) 



Taking into account condition (77), from non-negativity of p(z,s) it follows that qio > 0. By virtue of Eq. ( 26 1 , this 
condition implies, in turn, that q > and, consequently, due to Eq. (60 1, that go > 0. 



• If we choose e = ±1, then Eq. (77l implies that 6' < 0, while Eqs. ( 74 1 and (79 1 take the form of standard 
'dark solitons'. i.e. 



and 



p(z,s) 



F(0 = -/tanh 2 



9 °^ (s) tanh 2 



K\qo\ 



(80) 



(81) 



• For e ^ ±1, condition (77l and Eq. (79 1 take the form of standard 'grey solitons', in the following range of the 

- 1 < e < 1 and 6' Q < . (82) 



parameters e and 9' 0l i.e. 
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IV. THE EVOLUTION AND COLLAPSE OF THE CONTROLLED GPE STRUCTURES 



In this section, we carry out a stability analysis of our system, taking into account the control condition (59 I and the 
explicit dependence of the controlling parameter quj on the transverse scale lengths a x and a y through Eq. (26 1. For 



simplicity, we present the results for the simple case when the perpendicular solution is the product of only two single 
Hcrmite-Gauss modes, which implies for qio the expression (30 1 for the transverse mode (n, I). The generalization to 



the multiple modes' solutions may be somewhat lengthy, but is straightforward. 

We note that the characteristic spatial scales a x , <r y , and ?™^(s) are related with the coefficients of the restoring force 
u> x , ujy, and u) z , through four equations. These are the Ermakov-Pinney equations (23), the minimization condition 



for the control, Eq. (30 1, and the consistency condition of the amplitude and the phase of the parallel solution, Eq. 
PI, 



ds 2 



+ w 2 As)(ij (s) - 



Aa^sf 



= 0, 



J = x,y, 



(83) 



n.l / \ 

lib ( s ) 



$n Si 



<Jx (S) <Jy (S) ' 



(84) 



2 / \ n.l u 



.lib, 



(85) 



Obviously, if we wish to produce a controlled state in a BEC experiment, only two coefficients of the restoring 
force can be adopted arbitrarily, while the third is determined by the above constraints. Alternatively, we may 
adopt two arbitrary dependencies in the form Fi(u) x , lo v , uj Zi s) — 0, F2(u> x ,u) v ,u) z ,t) = 0, from which the quantities 
uij(s), j = x,y,z are uniquely determined with the use of the constraints ( S3 1-( 85 1 . The temporal evolution and the 
stability of the resulting coherent nonlinear state depends strongly on our choice of the temporal dependence of the 
confining potential Vt rap - Both stable and collapsing solutions can be obtained under different conditions, which is 
studied in more details below. 



Stable solutions 



First we study the important particular case when the perpendicular restoring force is stationary. For this case, we 
demonstrate the existence of stable (non-collapsing) nonlinear modes when the perpendicular solution contains only 
one Gauss-Hermite mode in each direction, see Eq. (21 1. 



The Ermakov-Pinney equation (23) is easily solved if ujj = constant, where j — x,y. The general solution for the 



perpendicular scale and for the phase functions, Eqs. (24 1 and (25 1, can be written as 



V^Vj.o + 4tJ M-0 + 8 7j,oWj<r4 sin^-s) + A;; - [4 ( 7 ? - uf) <rj + A;?] cos(2w jS ) 
°} ( s ) = ^ . (86) 



and 



7i (s) 



{8-f jfi uJ : jaj cos(2uj 3 s) + [4 ( 7 | - uf) a^ + Ag] sm{2u js)} 

4 ^VI,o + 4w M,o + 8 7i,o^i4o sin(2 WjS ) + A? - [4 ( 7 ? - % 2 ) + A 2 ] cob(2w 3 -s) : 



,k (s) = <f>j,k,0 



2fc+ 1 



tan 



27j,o4o , 4 7?o4o + A 



A,- 



tan (ujjs) 



tan 



An 



(87) 



where ct^q = crj(0), 7 j.o = 7/(0), and , 4>j,k,o — ^j,fc(0) are arbitrary initial conditions. 

We now restrict ourselves to a 2D geometry, d/dy = 0, where the analysis is particularly simple, but it can be 
easily generalized also to the 3D case. In a 2D case we have 



za (s) = 



1 



— i?n 



2*1 L^-i 



9o A c 



tan 



2 7^o< , 4 7 2 ct^ + a2 



A, 



tan (lo x s) 
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tan 



(89) 



and wc also conveniently rewrite the expression for a x , Eq. ( 86 1 , in the following form 



al (a) = DU 2 cos 2 [w x (s - S x )] 




(90) 



where D x and S x are two arbitrary constants. It is obvious that ct 2 (s) is non-negative, and that it may have periodic 
zeros only in the absence of spatial dispersion, A c = 0. Thus, the amplitude and the characteristic wavenumber of the 
solution, oc qijj oc l/a x (s) are regular functions, and the collapse does not occur. The corresponding variation of the 
coefficient of the parallel restoring force, w 2 , in a 2D regime with a single Hermite-Gauss mode, using Eq. (59 1, lo z is 
readily obtained as 



2 cos 2 [u x (s-S x )] -1 + W1 + 



X 2 . 



^ X D X 



cos 



i [2cj x (s ~ S x )} + V /1 + A 2 /(4^D4) + \J (2u x Dl 



{cos [2w x (s - S x )] + ^TTWW4Df)Y 

1 

cos [2u x (s - S x )] 



(91) 



The first term is positively definite, and the second is the sum of a harmonic function and a constant, of the form 
cos[2wj_(s — S±)] + C, where < £ < 1, which obviously features the change of sign. 




FIG. 1: Left: The temporal evolution ("breathing") of the perpendicular scale a x 

(dashed), found in a stable configuration, with \qo\ = 1,A C = 
.15, (fixko — 0, and for the lowest order Hermite-Gauss mode, 



characteristic frequency w z (s) in the 2D case, given by Eq. (911 
defocusing and vice versa. 



(solid line) and the parallel scale qio(s)/q 

1,5 = \ c \q \/^/2\§[\ = l,(j x = l,a x , = 0.607107, -y^o = 
k — 0. Right: The corresponding solution for the parallel 
There is a periodic change of sign, from focussing to 



The temporal evolution ("breathing") of the perpendicular scale c" 1 and the corresponding evolution of w 2 (s) in 
the 2D case are shown in Fig. [I] respectively. Hereafter, we fix for all figures: \qo\ = 1. 

The evolution of the 2D structures that in the parallel direction correspond to a bright, dark, and gray soliton, 
described by Eqs. ( [69] ) , ( |80[ ), and ( [79] ), respectively, is for the case of a ground Hermite-Gauss mode, k = displayed 
in Figs. 2|4 The solution features "breathing" (the oscillations of the amplitude and of the position) due to the 
oscillations of the perpendicular solution, as well as oscillations in the parallel direction, described by the quantity 
zq(s), Eq. (70l. The frequency of the latter oscillations is smaller than that of the "breathing", and they are not 

2l|4l 

=1) mode, is displayed in Figs. 



noticeable within the time span of Figs 

The evolution of the stable, first order (k 
ground mode. 



5|7 Its behavior is similar to that of the 
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t=o 



t=0.45 



t=0.9 



t=1.3S 







t=1.8 



t=2.24 



t=2.7 



t=3.14 







FIG. 2: The evolution of the lowest order (k = 0), 2D controlled, stable BEC state whose parallel component is a bright soliton. 
The initial position is calculated with Rq — 0.2, and the other parameters of the solution are the same as in Fig. [l] 



t=0 



t=0.45 



t=0.9 



[ 1.3? 







t 18 



t=2.24 



t=2.7 



t=3.14 







FIG. 3: The lowest order (k = 0), 2D controlled, stable BEC state whose parallel component is a dark soliton. The parameters 
of the solution are the same as in Fig. [T] 



B. Collapsing solutions 



Besides the stable solutions described in the preceding subsection, unstable controlled BEC states can also be 
generated, with a different choice of the confining potential Vt ra p- This can be demonstrated in the simple case when 
the parallel restoring force is adopted to be stationary, u z = constant, when the equation (59) can be solved as 



9id (s) 



1 



Ci COS (u z S + (fo) 



(92) 



where c\ and <^o are arbitrary constants. Then, in the 2D case {d/dy = 0) and using Eqs. (30 1, (23 1, (24), (|25j ), and 
(70 1, we can readily write down the parameters of the perpendicular, Eq. (22 1, and parallel, Eqs. (69 1, (|80h , and 
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1=0 



t=0.45 



t=0.9 



t=1.35 





t=l.B 



t=2.24 



t=2.7 



t=3.14 





FIG. 4: The lowest order (k = 0), 2D controlled, stable BEC state that in the parallel direction behaves as a gray soliton with 
e = 0.7. Other parameters are the same as in Fig. [l] 



f=0 



t=0.4? 



t=0.9 



t=1.3> 







t=1.8 



t=2.24 



t=2.7 



t=J.14 







FIG. 5: The evolution of the first order (k = 1), 2D controlled, stable BEC state whose parallel component is a bright soliton. 
The parameters of the solution are the same as in Figs. |2|4| 



(79 1, solutions, viz. 



P x (s) = [uj z cos (uj z s) + 7 Xj0 sin (uj z s)} , 
w 2 



Ix (s) = UJ, 



7x,o cos (u) z s) - uj z sin (u/ z s) 
7 Xj0 sin (ui z s) + lo z cos (u> z s) ' 



(93) 
(94) 



,k (s) — 4>x,k,0 — 



A c (l + 2fc) sin(w z s) 

4 al u z cos (u x s) + 7^,0 sin (u> z s) ' 



(95) 
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1=0.45 



1=0.9 



fcl.35 







t-1.8 



t=2.24 



t=2.7 



1=3.14 







FIG. 6: The first order (k = 1), 2D controlled, stable BEC state whose parallel component is a dark soliton. The parameters 
of the solution are the same as in Figs. |2| 1| 



M> 



t=0.4> 



1=0.9 



t 1.3? 







1=1.8 



t=2.24 



I 27 



t=3.14 







FIG. 7: The lowest order (k = 1), 2D controlled, stable BEC state that in the parallel direction behaves as a gray soliton with 
e = 0.7. Other parameters are the same as in Figs. |2| l| 



/ X 2 U! 2 

uj x (s) = uj z a 1+ tY [u g cos (u> z s) + 7^0 sin (u> z s)]~ 



qiD (s) = [w, cos (u; z s) + 7^0 sin (w z s)] 1 
(sj = i—r H sin (w z s) , 



(96) 
(97) 
(98) 

?ir> (s) <?oo-x,o^ 

where 6^ is defined in Eq. (29 1, while o x § = u x (Q), 7 x .o = 7r(0), 4>x,k,o = C£c,fc(0), and i?o are arbitrary initial 
conditions. 
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We note that the characteristic wavenumbers and the amplitudes of both the perpendicular and the parallel com- 
ponents of the BEC controlled structure are periodic functions of time, with the frequency co Zl which take infinite 
values twice during each period. The displacement from the equilibrium position, z (s) /<7id(s), oscillates with the 
same frequency, but remains limited. Such behavior corresponds to a cyclic collapse and recovery of the soliton-likc 
structure in both directions. However, the recovery does not occur in reality Our basic GPE is no longer valid for 
large amplitudes, such as those associated with a collapsing solution, and the new physical effects that arise under 
such circumstances are likely to perturb the system in a way that prevents its recurrence. 

This kind of collapse is well known for the bright solitons, Eq. (69), in the case of a focussing nonlinear term 

in the GPE, qio < and 8' > 0. Our analysis shows that the collapse occurs also for a defocusing nonlinearity, 
Qid > and 8' < 0, when the nonlinear solution is periodic in z. The wavelength of such train of nonlinear structures 
is proportional to 1 / ^/qm, and becomes infinitcsimally small in the singularities of qm, i-e. the train of coherent 
structures undergoes a collapse. 




FIG. 8: Left: The temporal evolution of the perpendicular scale a x 1 (solid line) and the parallel scale qm(s)/q (dashed), 

found in a collapsing configuration, with \q \ = 1, A c = 1,8 = A c |g |/y 2\8' Q \ = 1,lo x = l,a x ^ a = 0.607107, 7^0 = -15, 4> x .k,o = 0, 
and for the lowest order Hermite-Gauss mode, k — 0. Right: The corresponding solution for the perpendicular characteristic 
frequency u x (s) in the 2D case. Note that the collapse needs to be supported by a rapid (explosive) growth of the restoring 
force. 



The temporal evolution of the perpendicular and parallel scales, u" 1 and qn>(s), are displayed in Fig. [8] together 
with the coefficient of the perpendicular restoring force, oj x (s). The collapse of the 2D structures, that in the parallel 
direction correspond to a bright, dark, and gray soliton is displayed in Figs. [9 |ll| in the case of a ground Hermite- 
Gauss mode, k — 0. The solution features an explosive collapse within a limited period of time. The evolution of 
the unstable (collapsing), first order (k = 1) mode, is displayed in Figs. T2~fl4 Its behavior is similar to that of the 
ground mode. 

If there is an initial displacement in the parallel direction, i? 0j m t ne course of its collapse, the soliton-like 
structure will approach the equilibrium position z = 0. The collapse occurs roughly at the same rate in both directions, 
and that it needs to be supported by a rapid (explosive) growth of the restoring force. 

The collapse of the controlled BEC structures, that occurs in the regime ui z = constant and with an arbitrary 
dependence F 1 (ui x (,lu v , s) = 0, is essentially 2D. It is obvious from Eqs. ( 92 1 and (30) that the singularities of the 
effective parallel and perpendicular wavenumbers scale as qm(s) cx l/[a x (sja y (s)] cx l/(s — sc)- Thus, if we adopt 
Lu y = constant, we obtain a regular expression for a y , Eq. (86 1, i.e. the collapse would occur in the x,y plane, 
producing a thin ID string along y. Likewise, if we adopt the dependence a x — a y (which is a stronger requirement 
than lo x = u)y), the collapse in the perpendicular, x and y, directions will be the same, and much slower than in the 
parallel direction. In other words, the controlled BEC would collapse first into a pancake structure in the perpendicular 
plane, that would at a later stage and on a slower time scale, collapse also in the radial direction. 



V. CONCLUSIONS 



In this paper, we have used a newly developed analytical procedure to construct a wide class of localized solutions 
of the controlled 3D GPE governing the dynamics of BECs in the presence of a spatio-temporally varying external 
potential, where the latter is composed by a 3D parabolic time-varying potential trap plus a controlling potential 
to be determined self-consistently. This has been done on the basis of recent investigations [33J [33] . According to 
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1=0 t=0.21 t=0.42 t=U5 




9 9 9 



FIG. 9: The evolution of the lowest order (k = 0), 2D controlled, unstable BEC state whose parallel component is a bright 
soliton. The initial position is calculated with 7?o = 0.2, and the other parameters of the solution are the same as in Fig. [8] 



t=0 t=0.21 {=0.42 t=1.35 




9 9 9 



FIG. 10: The lowest order (k = 0), 2D controlled, unstable BEC state whose parallel component is a dark soliton. The 
parameters of the solution are the same as in Fig. [8] 

these results, we have found a class of solutions in the factorized form i/>(r, s) = 4>±( r ±i s)ip z (z, s), which allows us to 
decompose the 3D GPE into a pair of coupled partial differential equations (i.e. a 2D linear Schrodinger equation, 
governing the evolution of ip± and a ID NLSE, governing the evolution of ipz), plus a constrained variational conditions 
which, among all possible choices of V con t r , does not change the mean energy of the system and therefore minimizes 
the effects introduced by our control operation. We have extended our previous investigation [34 to a wider family 
of localized 3D solutions of the controlled GPE. In particular, we have found a wide class of localized solutions whose 
transverse profile is expressed in terms of breathing Hcrmitc-Gauss functions and the longitudinal one expressed 
in terms of breathing bright, dark or grey solitons. Furthermore, we have studied in details the properties of the 
controlled BEC states in the 2D case. It is demonstrated that they feature the breathing due to the oscillations of 
the perpendicular solution, as well as the oscillations in the parallel direction, arising from the initial displacement of 
the structure from the bottom of the potential well in the parallel direction. In the examples displayed in our Figs. 
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1=0.21 



t=0.42 



t=1.35 




t=1.8 



t=1.0? 



t=4.26 



t=1.47 




FIG. 11: The lowest order (A; = 0), 2D controlled, unstable BEC state that in the parallel direction behaves as a gray soliton 
with e = 0.7. Other parameters are the same as in Fig. [8] 
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t=1.35 







t=1.8 



t- 1.26 







FIG. 12: The evolution of the first order (fc = 1), 2D controlled, unstable BEC state whose parallel component is a bright 
soliton. The parameters of the solution are the same as in Figs. |9|ll| 



[T] [T4j the frequency of the latter oscillations was smaller than that of the " breathing" , and the parallel oscillations 
were not noticeable on the "breathing" time scale. 

The stability of the controlled structures, both bright and dark/gray, is governed by the shape and the temporal 
dependence of the potential trap. We note that the minimization condition for the control, Eq. (30), introduces an 



additional constraint for the possible experimental realization of the trap. The parallel and perpendicular coefficients 
of the restoring force are related, and they can not be both adopted arbitrarily. We have demonstrated that a stable 
solution can be obtained if the perpendicular restoring force is stationary, while the parallel force periodically changes 
the sign. In other words, the external parallel force needs to be switched periodically from the confining to dcconfining, 
in order to arrest the inherent collapse of the soliton due to the nonlinear interactions. Conversely, when the coefficient 
of the parallel restoring force is stationary, the controlled solution is unstable and undergoes a collapse within the 
finite period of time. To maintain the control throughout the collapse, the corresponding perpendicular restoring 
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t- t- 0.21 t=0.42 t=1.35 




9 9 9 



FIG. 13: The first order (k = 1), 2D controlled, unstable BEC state whose parallel component is a dark soliton. The parameters 
of the solution are the same as in Figs. |9|11| 

t=0 t=0.21 t=0.42 t-1 35 




FIG. 14: The lowest order (k = 1), 2D controlled, unstable BEC state that in the parallel direction behaves as a gray soliton 
with e = 0.7. Other parameters are the same as in Figs. [9 |ll| 

force needs to have a simultaneous singularity in time. In a 3D geometry, the collapse remains predominantly 2D. 
Thus, in the regime when one of the perpendicular coefficients is constant, the controlled BEC collapses into a ID 
string, while in a symmetric situation in the perpendicular plane, it collapses first into a pancake-like structure in the 
perpendicular plane, that will collapse also in the radial direction, at a later stage and on a slower time scale. 

We want to point out that this stability analysis is not exhaustive, but a more general stability analysis will be 
given in a future work. 
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